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THE USE OF MODERN TURBULENCE THEORY 
FOR CALCULATING EDDY DIFFUSIVITIES 


James C. Hill 

Laboratory for Theoretical Studies 


ABSTRACT 


The prospects are examined for obtaining analytical (as opposed to empiri- 
cal) solutions to turbulent transport problems. The direct interaction approxi- 
mation, a modern turbulence theory, is applied to the heat transfer boundary 
value problem and is seen to provide a reasonable generalization of the concept 
of eddy diffusivity. The calculation requirements are very great for present 
day computers. 
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SUMMARY 


The formal turbulent heat transfer problem is posed: given only the bound- 
ary data and certain statistical information about the turbulent velocity field, 
calculate the total rate of heat transport across the boundaries. The difficulty 
in calculating the joint velocity field-temperature field statistical moments 
stems from the closure problem of the dynamical moment equations. Of all the 
modem closure methods of turbulent theory, that one which seems most likely 
to have success in realistic heat transfer problems is the direct interaction ap- 
proximation (DIA) of Kraichnan. The DIA happens to be the exact solution for a 
model dynamical system and, consequently, cannot give the nonphysical behavior 
of other approximations not having positive definite probability density func- 
tionals. Furthermore, it properly allows the existence of an eddy diffusivity 
tensor in the case of a uniform mean temperature gradient. The latter can exist 
if the integral length scales of the turbulence are much smaller than any of the 
scales of statistical nonhomogeneity. 

The DIA equations involve the Green's function <G> for the mean tempera- 
ture field. In physical coordinates < G> is initially a delta-function and, conse- 
quently, the source of severe numerical difficulties. In wave-number coordi- 
nates these difficulties are avoided, as < G > is initially unity or at least diagonal 
with respect to wave-numbers in directions of nonhomogeneity. Unfortunately, 
the nonlinear convolution term requires six additional integrations for each 
direction of nonhomogeneity. This number can be reduced to two for a suitable 
expansion of the velocity covariance. 

In either coordinate representation the DIA places severe demands on the 
largest and fastest computers available. Because the DIA can be made ener- 
getically consistent for a finite mode system, it is recommended that the Fourier 
representation be used for computation in spite of the greater number of inte- 
grations compared to the physical coordinates representation. 
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THE USE OF MODERN TURBULENCE THEORY 
FOR CALCULATING EDDY DIFFUSIVITIES 
I. INTRODUCTION 

It is not possible at the present time to perform analytical calculations of 
turbulent transport problems with the same degree of confidence and rigor one 
would have for problems in ordinary molecular diffusion, heat conduction, or 
laminar convection. Unless the turbulent velocity fluctuations are known with 
absolute certainty, the mathematical problem may not even be solvable in prin- 
ciple. 

Consider the familiar heat conduction equation 

3 T 

R -kV*T = 0 


where T is the temperature and * is the constant molecular thermal diffusivity 
(or thermometric conductivity). Many analytic solutions are known for this equa- 
tion. The situation is less hopeful with the addition of a laminar convection 
velocity l (x, t) 

3 T 

T- + U- 7T-xV 2 T = 0. 

o t 


Only a few solutions for special IJ are known, although approximate solutions 
may be obtained numerically for a wide class of problems. When li takes on 
random or turbulent values u, or LI = <L> + u, where the angular brackets de- 
note the average or known part, the situation is almost hopeless, as the random 
coefficient u in the equation 

3 T 

+ <l)> • TT * x V 2 T = - u • 7 T (1) 


induces a random part 9 of T, or T = <T> + 9. 

The reason that u is considered unknown and not merely complicated lies 
in our inability to solve the Navier-Stokes equations with sufficient accuracy. 
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Because of their nonlinearity and nonlocalness the equations of motion are 
extremely sensitive to initial conditions. An infinitessimal error will generate 
at some time a completely different flow. This is related to the observation that 
fluid particle trajectorie*' *nd statistically to wander apart. Alternatively, we 
can look at this from .int of view of Heisenberg; 1 a fluid is a mechanical 
system with an infinite number of degrees of freedom, all of which must be speci- 
fied initially to determine the motion uniquely. Only when the viscosity is high 
enough (or the Reynolds number low enough) will the number of degrees of free- 
dom reduce to one (laminar flow), for which L can be calculated unambiguously. 

Although an exact solution to the turbulent heat transfer problem is hopeless, 
a statistical one may not be. Here we will be concerned with the statistical 
theory of turbulent convection, in which statistical averages, denoted by angular 
brackets, are taken over an infinite collection or statistical ensemble of flows. 
Each member flow[n] is called a realization. We have, then, 



for the mean of U and similar definitions for averages of other quantities. 

Typically, eddy diffusivities and eddy viscosities do not arise in modem 
theories of turbulence and turbulent convection. They were introduced by 
Boussinesq 2 in 1877 and given physical basis by Prandtl 3 with the crude concept 
of mixing length. These approaches have been more successful so far than any 
modem analytical theories because of the presence of adjustable empirical con- 
stants, rather than their being good representations of the physics (Reference 4, 
p. 275). 


H. THE PROBLEM 

To define the problem let us consider as an example the case jf heat con- 
vection by an incompressible single-phase fluid in a region $ having boundaries 
®. We assume constant physical properties with no viscous dissipation. The 
problem is to determine the heat transfer across ft, <q>, given only the initial 
and boundary temperature distributions T 0 (x) and T B (x, t), the Prandtl number, 
and the Reynolds number. We assume that the Reynolds number determines all 

the statistics of l) and that these are known. Because <q> = - k 3 — < T> n , the 

dn B 

problem reduces to finding the average temperature distribution <T(x, t)> 
throughout ft. 
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The fundamental difficulty of solving the statistical problem is illustrated by 
the moment formulation cf the problem. Rewrite Equation (1) as 


£ T = - u • V T 


where £ is the nonrandom operator 


£ * 



7 - 


K 



Averaging Equation (1), we get 

£ <T> -■ - 7 • <u 6>> (2) 

which involves an undetermined quantity, the velocity-temperature correlation. 
This equation was given by Kampd de F6riet s in 1937. A closed equation for < T > 
may be formed from Equation (2) with the introduction of an "eddy diffusivity" e 
defined by* 

<u d> = - e 7 <T> 


and which must be evaluated empirically. This approach is very crude, inasmuch 
as Equation (2) must be solved— empirically, at least— in order to determine e . 
Alternatively, we may use Equation (1) to write eq”ations for the higher order 
undetermined moments that appear in Equation (2) and similar equations. An 
unclosed hierarchy of moment equations results, as in writing an equation for 


*ln general £ must 


be considered a second order tensor (Cf. Equation (6) and Reference 4, p. 279). 
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the higher order moment that occurs in the preceding equation, the nonlinearity 
of dynamical quantities always produces another unknown. We can obtain in the 
present problem the hierarchy 

£ <T> = - 7 • < u ti> 


£ < u ' 6> - - V • [ '" u u' 6’>t<uu'><T>] 


£ < u' . . . u" 0 > 

v _ ; 

T 

n 


- 7 


II u 


u" d> 


r+T 


< u ' . . . u"> < ii U > + 

' V ' 



n + 1 



where we have used nonsimultaneous and noncoincident u's (denoted with primes), 
that are transparent to the £ and 7 operators, in order to avoid the Navier- 
Stokes equations. In this formulation certain space-time arguments are made to 
coincide before substitution into the lower order moment equation. The task of 
suitably truncating this set of equations to a deterministic set is called the 
closure problem. 

The closure problem of the moment equations illustrates the fundamental 
difficulty of turbulence and turbulent convection theory. In each realization T is 
a functional of the infinitely complicated function U (x, t). The infinite complexity 
of the turbulence t (x, t) and its effect on T (* , t) cannot be removed simply by 
averaging. 6 This difficulty persists regardless of the formulation of the problem. 
In the moment formulation, besides the closure problem of the hierarchy of mo- 
ment equations, iteration series (expansion in terms of a turbulent P6cl6t num- 
ber) do not converge. In the probability functional or characteristic functional 
approaches, multimode interactions are inseparable in wavenumber coordinates 
and nonlocal interactions are nondiagonalizable in physical coordinates. 


fil. APPLICATION OF A MODERN TURBULENCE THEORY 
A. Status of the Theory 

How will modern turbulence theory help us? Usually the same approaches 
are used for turbulent convection as for describing turbulence itself, although 
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one must be careful, because the nature of the nonlinearities are different. Most 
of these approaches are not satisfactory. This is not surprising, since almost 
all of them are based in some sense on perturbation expansions about the case of 
zero turbulence. Any success with these in real problems must be considered 
luck.* 

One lucky theory is the direct interaction approximation (DIA). It is probably 
the most satisfactory of modern closure procedures as far as its ability to give 
sensible results is concerned and is the one we will discuss here. It was invented 
by Kraichnan 8 9 in 1958, given for the initial-value scalar convection problem in 
1961, 10, 11 and can be easily generalized to the arbitrary boundary- value heat 
transfer problem.* It has been derived in several different ways. 

Although there is some disagreement over the extent to which the DIA is a 
"rational approximation", 14 * 1, we believe it to be a good approximation, or at 
least, one of ihe best available, because of its ability to satisfy a number of con- 
sistency requirements and properties of the true solution. These desiderata for 
closure theories should be satisfied by any turbulent convection theory that is 
c pectcd to give go>>: results. The following list contains some of the more im- 
portant reqvirements.* Unfortunately, no existing closure scheme satisfies them 
all. 


(1). Global Conservation Laws. The DIA and the true system g.ve the same 
result when averaged over B. The DIA preserves global conservation of thermal 
energy 12 

-77 I <T> dB * f <T> <li> • d» - « I ▼ <T> • d:B 

dt Jfl Jjg 


and of an "entropy" 


I <T 2 > dB + f <T 2 ><U> • dB = k f V <T 2 > • d8 - 2k f < |VT| 2 > dB. 
dt J B J B h J % 


*This pessimistic view is due to Kraichnan. 7 Some of the approaches are referenced in Appendix 
A. 

f Kraichnan 12 had considered the boundary -value problem but did not properly take into account 
Equation (3). See also Reference 13. 

$ This list is taken from Kraichnan 7 and from Orszag. 17 
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In these equations we have assumed that fluctuations in u disappear on®. 


(2). Realizability. All probabilities associated with the random variables T 
and l must be positive. This requires positive-definiteness of powder spectra, 
oounds on the skewness, and inequalities for all the other moments. 1819 The 
DIA satisfies this requirement because it is the exact solution to a model system 
in which the infinite ensemble of realizations is broken up into an ensemble of 
collections of realizations where u and T randomly interact. This situation is 
analogous to dime-flipping experiments in w'hich the outcome of one toss depends 
on all the others. The behavior of T in realization [n] is given by the equation 12 

£ j[ n l T _ V 1 <J)n „f k l . y j[ m ] 

/ k, m 


where 




k . m 


N 


“2 


r-.y 


exp * 


N 


[p(k - n) f >• (m - n)j 


“"p+r. 


N is the number of realizations in a collection, and the $ 's are restricted by 

k! - i 


= i 0- 


or y 


0) 


■p, P 


but otherwise take on random phases throughout the collections. 

(3). Galilean Covariance. The DI equations do not transform properly under 
Galilean coordinate transformations (convection by a uniform velocity) in homo- 
geneous turbulence. 20 This is perhaps the mair. failing of the DIA. The invari- 
ance of the dynamical equations to uniform convection should be retained to 
properly describe the convection of small scales by large scales, as necessary 
for the derivation of Kolmogorov's law* in ordinary turbulence theory. 

f 

| 

£ 
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(4) . Stochastic Irreversibility . It is important that relaxation effects be 
described correctly, at least qualitatively. In particular, we believe that the 
statistical equations should exhibit the property of a mechanical system with a 
huge number of coupled degrees of freedom to approach an asymptotic statistical 
state that is independent of the initial conditions. 17 The relaxation is controlled 
by two factors: one is the explicit irreversibility of the conduction term which 
provides a sink for the dissipation of temperature fluctuations and which can be 
diagonalized or localized by using wavenumber coordinates. The other is the 
nonlocalness in both physical and wavenumber coordinates of the convection term, 
which, we believe, tends to drive the system irreversibly towards equip# ”tition 

of the temperature fluctuations (destruction of spatial correlations).* The DIA 
retains both mechanisms. 

(5) . Quantitative Contact with Dynamical Equations. Besides the qualitative 
features we have mentioned, it is important that the DIA be in some sense an 
approximation to the dynamical equations (Equation (1) plus the Navier-Stokes 
equations). This condition is met, as the DIA becomes exact as a perturbation 
solution in the limit that the turbulence Lecomes weak — that is, in the limit that 
the right hand side of Equation (1) may be considered a perturbation. 10 Require- 
ment (1) also serves to provide some contact. 

B. The Direct Interaction Equations 

Before writing down the DI equations, we need to define the Green's function 
for the temperature field. The Green's function G is defined to be the solution 
of Equation (1) for the special case of delta function initial conditions and zero 
boundary conditions — that is, G is the response to a unit pulse of temperature 
when the initial temperature is otherwise zero and the boundaries are infinitely 
conducting. Hence we have two problems, which we write symbolically as 

G - problem 


Jig - o 


Gb = 0 

itial ~ ^ ( ) 


T - problem 

fR: % i = 0 

S: T = T b 

10: T = T 0 


* 


These ideas have recently 


been incorporated in a phenomenological approach by Leith, 


21 
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where 


Jl 



L • V - k V 2 


The temperature field T can then be recovered by integrating over the initial and 
boundary data weighted by the Green's function, 22 which result we write sym- 
bolically* 


T = f T. G dJt - « f f T, »G • d:8 dt . 

k 0 Jo J* B 


( 3 ) 


We actually need <G>, however, for which there is an almost identical closure 
problem as for <T>. To simplify matters we use a maximal randomness condi- 
tion that T 0 and T B be statistically independent of the u field to all orders of u. 13 
This guarantees that in Equation (3) G will be independent of T 0 and T B , and 
hence <G> can be determined independently of <T>. 

The DIA yields the following closed integrodifferential equation for <G> in- 
volving velocity field statistics to only second order 10 

£ <G> = [▼ ■ <uu'> <G>] * [• V <G>] . (4) 


Here the asterisk denotes i space-time convolution integral. Equation (4) may 
be integrated over T Q and T B (Equation 3) to obtain an equation for <T>, which, 
by comparison with Equation (2), yields 12 

<u6> = - [<uu’> <G>] * [ • V <T> ] . (5) 


This goes to the form 


<u e> = — c • v <t> 


( 6 ) 


In Equation (3) the del operates on the source point. Equations (3)— (7) are written out more 
explicitly in Appendix B. See also Reference 13. 
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in the special case that 7<T> is approximately constant, where t is effectively 
an eddy diffusivity tensor Ut 12 


E 



<uu'> <G> d!Rdt . 


( 7 ) 


We see that the DIA provides a generalization of the classical ideas of eddy 
diffusivity that is consistent with the physics of the problem, that turbulence is a 
somewhat nonlocal phenomenon. In general the turbulent temperature flux <us> 
is not directly proportional to V<T> but depends on the entire temperature field 
and on the turbulence throughout An eddy diffusivity does exist only if V<T> 
is uniform over the region in space-time in which the velocity is correlated. This 
can occur at steady state in situations where the integral length scale of the tur- 
bulence is much smaller than any scales of statistical nonhomogeneity or, crudely 
speaking, if the turbulent eddies are much smaller than the dimensions of the 
flow system. The eddy diffusivity is tensorial and goes to the scalar form only 
if the turbulence is isotropic. These restrictions and generalizations on the 
existence of a classical eddy diffusivity seem to be qualitatively correct, and are 
analogous to the conditions for existence of a molecular thermal conductivity in 
the statistical theory of transport in gases. 

At this point it is natural to ask what sort of physical picture is represented 
by the DIA? The answer is difficult, as from what we have said the DIA seems 
to be an arbitrary mathematical model involving unrealistically coupled realiza- 
tions. However, consider the case of homogeneous turbulence for which we write 
the Fourier transform of Equation (1) as* 


T(k) + Kk 2 T(k) = - i k 


L 

k=k‘+k" 


G(k') T(k"). 


( 8 ) 


The model equation for which the DIA is the exact statistical solution is 
T(k) + Kk 2 T(k) = ik- 4 > (k, k', k”) u(k') T(k') 

k = k'+k " 


'‘This form is suitable if we let the T and u fields be periodic over arbitrarily large box-volumes. 
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where M = 1 , <£(any argument = 0) = 1 , <?(k, k',k") = <£*(k", -k', k) = 

<^3 (k» k' » k') = <P* (- kt ~k'* “k") and otherwise the phases of <p are random. The 
model, then, replaces the true wavenumber interactions by fictitious interactions 
whose strengths remain identical to those of the true system (because of the 
constraints on <p) but whose phases are randomized. The DIA gets its name by 
treating each elementary triad wavenumber interaction (a direct interaction) in 
the sum term of Equation (8) as a perturbation against the background of all the 
indirect interactions, and assuming that only the sum of the DI's survives the 
summation in Equation (8). With this mechanism excitations are thus relaxed by 
dynamical interaction with the full turbulent field, and not just by conductive 
decay. 


IV. COMPUTATIONAL CONSIDERATIONS 

We must be able to solve Equation (4) for <G> to obtain any practical infor- 
mation from the DIA. In order to test the approximation we will consider a brute 
force calculation before attempting such ad hoc simplifications as lower dimen- 
sional models or special assumptions on the form of <G>. 

The usual procedure for numerically approximating an integrodifferential 
equation is to set up finite difference equations on a space-time or wavenumber- 
time mesh and to march the solution forward from the initial time layer. These 
difference equations should be consistent with Equation (4) as the mesh density 
increases and should preserve the desiderata that we listed previously. We usu- 
ally have to resort to a posteriori evidence of numerical stability and assume 
convergence to the true solution as a matter of faith. 

The physical demands presented by Equation (4) on existing high-speed elec- 
tronic computers are nearly excessive. Huge amounts of storage are needed be- 
cause <G> has up to 8 independent variables. Great speed is necessary for 
evaluation of the convolution integral term, which requires a 4-fold integration 
(or summation) in x-space or up to a 22-fold integration in k-space. 

We have found severe numerical difficulties in attempting to calculate <G> 
in physical coordinates (x-space). 13 Most of these are associated with the fact 
that <G> is initially a delta function and is difficult to approximate numerically. 
The symptoms are large oscillations that decay no faster than <G> but which 
can be reduced by decreasing the time increments. The oscillations are induced 
at small times by truncation of the differential operators in Equation (4) and can 
be important at larger times because of the time convolution. An equally serious 
difficulty is failure to preserve the consistency requirements listed earlier. 
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Solution of Equation (4) in wavenumber coordinates (k-space) is much more 
attractive. First, the same set of consistency requirements can be satisfied — 
even with a finite mode system. Second, the initial singularity of <G> is avoided, 
because its initial value is unity for all wavenumbers or at least diagonal with 
respect to wavenumbers in directions of nonhomogeneity. Unfortunately, more 
integrations are required for the convolution term in nonhomogeneous turbulence. 
In going to k-space, for each dimension of nonhomogeneity each gradient operator 
produces one additional integration and the < uu‘> <G> product produces 4. The 
total number (seven) can be reduced to 3 for a suitable Fourier expansion 
of <uu‘ >. 

At the present time we are completing the calculation of an eddy diffusivity 
for the case of a steady isotropic turbulence and a uniform mean temperature 
gradient.* In this calculation <G> is a function of 2 independent variables, and 
the convolution requires 3 integrations (5 nested FORTRAN DO- loops). For com- 
parison there is one data point for the measurement of velocity-temperature 
correlation in a grid-generated turbulence heated to a self -maintained uniform 
temperature gradient. 34 The space-time velocity correlation coefficient may be 
estimated from the data of Favre 25 or Frenkiel and Klebanoff. 26 

The next problem to consider is steady heat transfer to fully-developed 
turbulent channel flow, for which <G> is a function of 5 independent variables 
and the convolution requires 6 integrations (a reduction from 10) — 11 nested 
DO-loops. Once <G> is obtained, the Nusselt number for a variety of boundary 
conditions may be determined. One case, the flow between parallel plates at 
different temperatures, may be compared to the data of Page, et al. 27-29 Another 
case is the turbulent Graetz problem. We expect the channel flow problem to be 
very difficult. 


V. CLOSING 

Although progress is^ being made towards the goal of solving turbulent 
transport problems analytically, the picture looks somewhat grim at the present 
time. Perhaps the situation will improve with the next generation of high-speed 
computers or when simplified approximations to the DIA are found. 


*This problem was invented by Corrsin 23 in 1952. 
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APPENDIX A 


Recently, there have appeared several reviews and comparisons of some of 
the analytical turbulence and turbulent convection theories (References 6, 7, 16, 
17, 18 Ch. 3, 30, 31), Since no really comprehensive and critical survey has 
appeared, closure techniques are still in the hands of specialists (except for the 
mixing length theories). Here we do not attempt a complete list of techniques 
but only indicate the diversity of approaches — many of which have not been 
carried out to completion, and some of which are merely alternate ways of 
setting up the same problem. 

We regard the more popular approaches originating before 1958 as the 
classical ones. These include — besides the phenomenological methods 4,32,33 — 
the moment discards, 34,35,36 * various types of quasinormal or fourth cumulant 
discard, 32,37 * 44 and the maximum dissipation hypothesis of Malkus. 45,46 

The level of sophistication suddenly increased with the inception of 
Kraichnan's direct interaction approximation. 8,9 Kraichnan also developed a 
random coupling model 10 for the DIA, an ad hoc modification of a Lagrangian 
form of the DIA, 47 and some variations. A great number of papers have been 
generated from the DIA (the present one included) — too many to list here. 

Work has continued on the use of characteristic functionals and functional 
probability formulations (Hopf, 48 Lewis and Kraichnan, 49 Rosen, 50,51 Edwards' 
random phase approximation, 52 Herring's self-consistent-field approximation 53,54 ) 
and the setting down of BBGKY-type hierarchies. 55 * 57 

Recently, Wiener- Hermite expansions of the random variable have been 
examined. 58 * 64 In these investigations the expansions were made about "white 
noise", although expansion about "red noise" would seem preferable for 
describing relaxation effects. 

There are still other approaches. Soon after the DIA appeared, Wyld 65 and 
L. Lee 66 derived consolidated perturbation expansions for the triple correlations 
in diagram notation. (See also Kraichnan 10 and Edwards. 52 ) J. Lee 31 has modi- 
fied the quasinormal approximation in the scalar mixing problem by comparison 
with the DIA. Orszag 17,18 has used "inequality preserving" closures, based on 
consistency requirement #2, for the cumulant equations in homogeneous turbu- 
lence. Kamp6 de F6riet 67 has set up the Gram-Char lier approximation for 
bivariate distributions, intended for application to homogeneous turbulence. 


* 


There have been several applications to shear turbulence, heat transfer, etc. 
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This short list does not exhaust the number of possible approaches. For 
example, variational and information- theoretic methods that have found use in 
other fields of statistical physics (e.g., Jaynes 68 ) have only been discussed 
sparingly. 8, 17 69 The possibility of singular perturbation expansions about 
infinite Reynolds number should also be considered — requiring first a theory 
for infinite Reynolds number turbulence involving, say, distinct separation of 
energy and dissipation spectra. Finally, because of dynamic passivity of T and 
the absence of a feedback effect on errors in u, the use of model random fields 
u should be investigated for "exact” solution of Equation (1). 


13 


APPENDIX B 


The equations in the text that involved G were only written symbolically. 
More accurately, let G(x, t | x', t') be the temperature at the point x e JR at time t 
for the response to a unit pulse of temperature at the point x' e 3? at the time 
t' < t . The equation for G is 

Jt(x, t) G(x, t | x', t* ) = 0 for x, x' e iR and t > t ' 

with the initial condition 

G(x, t* | x\ t') = S(x - x') for x, x' eiR 

and boundary condition 

G(x, t | x', t') = 0 for x or x 1 e$. 

The last condition is sufficient for reciprocity. Equation (3) is then 


T(x, t) - f d 3 x' T (x') G(x, t( x', 0) - k f dt' f d 2 x' T B (x', t') 3-7 G(x, t|x', V ) 

J iR Jo n 

where n is the outward unit normal on % . Equation (4) becomes 

£(x,t)<G(x, t|x',t' )> = {* dt" fd 3 x"7 * <uu"> <G(x,t|x",t")> • V x „<G(x",t"!x‘, t')>, 

Jf 


Equation (5) is 


u(x, t) <9(x, t) > - - f dt' f d 3 x' < uu'> <G(x, 1 1 x', t') > • V , <T( x'. t' )> , 

Jo 


and Equation (7) is 


*(x, t) 


f dt' f d 3 x' <uu’> < G(x, 1 1 x', t') > . 

Jo h 
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